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We propose a metapopulation version of the Schelling model where two kinds of agents relocate 
themselves, with unconstrained destination, if their local fitness is lower than a tolerance threshold. 
We show that, for small values of the latter, the population redistributes highly heterogeneously 
among the available places. The system thus stabilizes on these heterogeneous skylines after a long 
quasi-stationary transient period, during which the population remains in a well mixed phase. 
Varying the tolerance passing from large to small values, we identify three possible global regimes: 
microscopic clusters with local coexistence of both kinds of agents, macroscopic clusters with local 
coexistence (soft segregation), macroscopic clusters with local segregation but homogeneous densities 
(hard segregation). The model is studied numerically and complemented with an analytical study 
in the limit of extremely large node capacity. 


I. INTRODUCTION 

Modern societies are often faced to segregation; dictated by race, religion, social status or incomes differences, 
it represents a major issue whose outcome can range from social unrest, to riots and possibly to civil wars. The 
understanding of the rise of such phenomenon has thus attracted a lot of attention from economists, politicians and 
sociologists [DU- 

In a couple of papers written in the late 60s, Thomas Schelling I3IS] proposed a stylized model to describe the 
unset of segregation and pointed out a counterintuitive but widely observed result: a well integrated society can 
evolve into a rather segregated one even if at individual level nobody strictly prefers this final outcome. Even when 
individuals are quite tolerant to neighbors of their opposite kind, allowing them to relocate themselves to satisfy 
their preferences - namely maximize their perceived local fitness/utility - will make segregation to emerge as a global 
aggregated phenomenon not directly foreseen from the individual choices. 

Since the pioneering works of Schelling the model has attracted the attention of the community of physicists and 
mathematicians, interested in its simplicity and in the emergent behaviors recalling models such as Ising and Potts 
ones IHHIS]. In its simplest form the model proposed by Schelling considers a population composed by two kinds of 
agents sitting on the top of a regular lattice (or a 1 dimensional ring), each site being able to receive at most one agent. 
The latter being allowed to hop to a new empty lattice site once unhappy, that is once the fraction of agents of her 
opposite kind in her Moore neighborhood, is larger than a given tolerance threshold. The striking result by Schelling 
is that segregation will emerge for tolerances slightly larger than 1/3, well below the “more natural bound” 1/2. 

Several variations have been performed onto the original Schelling model (see [7] for a review), for instance different 
ways of computing the agent utility [15] and rules according which agents do move (swap agents positions, restrict 
moves to increasing fitness locations [3 HO], to fixed distance places or following the path of an underlying network ED- 
All such models contains essentially two parameters, the tolerance threshold and the fraction of empty spaces; the 
latter having received much more attention because it is responsible for a phase transition in the model laiiniiisHis]. 

Recently the Schelling model has been improved in realism by considering a metapopulation scheme daniiii]. 
Nodes have a large but finite carrying capacity representing thus a bunch of residences in a district whose spatial 
extension can be neglected (well mixed population); happiness and thus willingness to move is thus conditioned by the 
fraction of agents of the opposite kind inside a given node, with respect to the nodes occupancy (1-node fitness) [T6] . 


II. THE MODEL 

In this paper we elaborate further in this direction by considering a metapopulation version of the Schelling model 
where two kinds of agents, say Red and Blue, can move across N nodes, each of which can receive at most L agents ^ 
(carrying capacity) at any given time. We assume there are in the model pLN vacancies, i.e. empty spaces, and an 
equal number of Red and Blue agents. 


^ To make things simple we assumed a constant carrying capacity for each node, but of course one can improve the mode by considering 
a different value of Li for each node i. 
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The considered model resembles to the one proposed in m, the main difference being that in our case agents don’t 
have any information about the selected destination, therefore also moves that decrease or leave invariant the fitness 
are allowed, we define such case weak liquid version of the Schelling model, the liquid case referring in the literature 
to agents’ moves for which the fitness does not decrease [9]. A second difference is that in [16] the agent fitness is 
computed using single node informations, namely the number of agents in a given sites, on the contrary we reintroduce 
as Schelling originally did, the concept of spatial proximity [6], agent fitness takes into account the number of agents 
in a given node and in the neighboring ones, namely nodes at distance 1 from the current node (local fitness). 

Nodes are assumed to be arranged in regular lattices, as initially assumed by Schelling, but the model can be easily 
extended to complex networks. The local update rules is defined as follows. The neighborhood of an agent is given by 
the topological neighborhood of the node where she lives, including the latter. At each time step an agent is selected 
and her fitness is computed as the fraction of agents of the opposite kind of her, living in her neighborhood with 
respect to the total number of agents living in the same neighborhood. Mathematically, assuming she is a Blue agent 
living in node i, then her fitness is given by: 


/f = 


E 


jei 




( 1 ) 


where X = is the number of agent of X-kind in node j, and we used the notation j G i to denote all nodes 

j belonging to a neighborhood of node i, including the latter, that is the set of nodes at distance smaller or equal to 
1 from i. 

As in the original Schelling model, agents are unhappy if their fitness is larger than a given tolerance threshold^ 
e G (0,1), hereby assumed to be the same for all agents of both kinds; to reduce their uneasiness agents move by 
choosing uniformly at random another node k not completely full, i.e. < L: 

if/f > e ^ agent A leaves node i. (2) 


The case for a Red agent is similar. 

One time step is the random selection with reinsertion of (1 — p)NL agents. We define the convergence time, to be 
the time needed for the system to reach the equilibrium, namely once no agents will move anymore. 


III. RESULTS 

We hereby present the numerical analysis of the proposed model once the underlying network is a regular lattice 
with periodic boundary conditions and each node has 4 neighboring nodes. The system is initialized with pNL 
vacancies, (1 — p)NL/2 Red agents and the same number of Blue ones, uniformly random distributed among the 
N = 400 nodes. The carrying capacity has been fixed to L = 100 and we check that the initial conditions satisfy the 
local constraint nf^nf<L for all i. Throughout the paper the emptiness has been fixed to p = 0.9. 


A. Single node properties 

The aim of this section is to present the local properties of the system, namely at the level of single nodes. The 
metric we used is the average value, over all the nodes, of the node magnetization: 

= -E 

N ^ 

i 

small values of (/i) mean that, on average, each node is populated by the same number of agents of both kinds, while 
large values are associated to nodes filled with agents of only one kind. 

For e < 0.5 we observe (see Fig. [^panels A and B upper plot) that asymptotically (p) ^ 1, meaning that the 
system stabilizes into a frozen state where local segregation is present in all nodes: each node contains only agents 
of one kind. Let us observe (see Fig. [^B upper plot) that the same behavior is also present in the simplified model 
where the fitness is calculated on the single node, as done in m, and thus it is intrinsic to the displacement dynamics 
and not to the way the agent fitness is computed. We also notice (see Fig. [^panels A and Fig. 6 that as e decreases 
toward zero, the time needed to reach the frozen state gets longer, going to infinity in the limit e ^ 0 . The system 
exhibits thus a transient phase where it remains stuck for very long time into a quasi-stationary non-segregated state 
(see red circles and orange stars curves in Fig. UK)- 


nt 


nf + nf 


(3) 




3 



FIG. 1: Single node behavior. Panel A: average magnetization {/j) = ^ ^ as a function of time for a single generic 

simulations. For e < 0.5 the magnetization goes to 1 indicating the total segregation inside each node. Panel B: asymptotic 
average magnetization (upper plot), average maximal node population maxi(n^ + nf) (middle plot) and fraction of empty 
nodes fempty (lower plot) as a function of e. For e < 0.5 the asymptotic average magnetization is constantly equal to 1, showing 
the robustness of the local hard segregation against replicas. As e ^ 0 agents of the same kind tend to select and fill all the 
compatible - with respect to the happiness - nodes, creating thus an important fraction of empty nodes. The black points 
correspond to the model with local fitness (namely distance-1 nodes) while the grey ones to the simplified model with 1-node 
fitness (namely distance 0 nodes). As we can see the magnetization behavior doesn’t depend on the way the fitness is computed. 
Panel C: distribution of nodes total population (n = tt-a + tib) for few values of e = 0.28, 0.3, 0.35, 0.5, 0.7. Averages shown in 
Panels B and C have been obtained performing 100 replicas of the model. 


A second fundamental self-organized phenomenon emerges for e < 0.5, the initial homogeneously distributed popu¬ 
lation organizes itself into an heterogeneous state across the network nodes (Fig. panel C lower plots): most of the 
nodes contain ^ 10 agents, while very few nodes have as much as ^ 100 agents, recall that L = 100 is the maximum 
node capacity. 

Observe that such asymptotic distribution is correlated with the time the system spends in the quasi-stationary 
non-segregated state, the longer this time the more the population distribution across nodes moves from a Poissonian 
distribution (Fig Impanel C upper plots) to a power law (Fig. panel C lower plots). The maximal node population 
i^max = uf) iucreascs for e ^ 0 (Fig. [l panel B middle plot). Notice that the maximal node capacity 

(L = 100) is never reached. At the same time we observe the formation of a relevant fraction of completely empty 
nodes, i.e. nodes for which nf + nf = 0 (Fig. panel B lower plot). 

The transition to the local segregation at e = 0.5 is an important issue of the metapopulation version of the Schelling 
model that can be easily explained once the fitness depends only on one node, with the following argument. A mixture 
of agents of both kinds in each node, determines a global happy configuration if and only if: 


Vi 


n\ 


n\ + n\ 


< e 


n\ + n\ 


< e. 


( 4 ) 


but it is straightforward to observe that this equation admits a solution, for which > 0 and > 0, only for 
e > 1/2. Hence for e < 1/2 the above relations cannot be satisfied and thus unhappy agents start to move, increasing 
the happiness of agents of the same kind and decreasing the happiness of agents of opposite kind. The net result of 
such behavior is segregation. 


B. Global properties 

To analyze the spatial structures, thus beyond the single node, we define a cluster based on node’s majority, more 
precisely two linked nodes belong to the same cluster if they both are characterized by the majority of agents of the 
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same kind 


i,jeC if {nf > nf and nf > nf} or {nf < nf and nf < nf} . (5) 

If a node share the same (positive) number of Red and Blue agents, then it will be considered part of the interface. 
An edge between two neighboring nodes i and j is considered interface if {nf — nf){nf — nf) < 0. Hence a cluster 
is made by nodes while an interface can contains both nodes and edges among them. 

This indicator allows us to show that (Fig. panels B and C) for e < 0.6 macroscopic percolating clusters are 
always formed, the average size of the largest cluster oscillating between O.AN and O.bN and at the same time the 
average sizes of the first and second largest cluster added together cover more than 75% of the available nodes. This 
critical threshold is the same observed in the original Schelling model. Notice that for e = 0.5 the first and second 
cluster cover almost completely the lattice, containing more than 95% of nodes. From the behavior of (Smax) foi* 
0.5 < e < 0.6 and the results of Fig. [^Panel B, one can conclude that the largest clusters contain a majority of agents 
of the same kind, but different agents can coexist in the same node. We call this scenario soft segregation because it 
allows a small mixing in the population. For e < 0.5 each spatial cluster contains only one kind of agent, resulting in 
a hard segregation of the population. For e ^ 0 clusters of empty cells are created increasing the distance between the 
two monochromatic structures by interposing interfaces (Fig. [^panels B lower plot and white nodes in Fig. [^panels 
A3, A4 and A5). 

Looking at the temporal growth of the interfaces (Fig. [^panels D), we can observe another remarkable self-organized 
phenomenon: for 0.3 < e < 0.6 the size of the interface, monotonically decreases in time as t~^^^ ( 2 : = 4 for e = 0.6 
and 2 : = 3 for e > 0.5). This is the typical signature of a coarsening phenomenon, that has already been observed in 
the classical Schelling model HO]. On the other side, for low values of the tolerance where the system remains for long 
time in the quasi-stationary state (for instance e < 0.3), the size of the interface has a slowly increasing phase (red 
curve Fig. Impanel D), corresponding to the formation of temporal segregated domains followed by an abrupt decrease 
when the system reaches the local segregation equilibrium. In this case the domain formation mechanism cannot be 
ascribed to the coarsening framework: the spatial segregation indeed is reached through a mechanism of aggregation 
of agents around the high density instabilities that allow the system to exit the quasi-stationary state. Once the 
clusters are formed, the complete equilibrium is reached as the interface between Blue and Red zones becomes empty 
(See Fig. 7 for a detailed explanation). 

While the creation of heterogeneity distribution of agents among nodes is due to the metapopulation mechanism, 
the formation of the monochromatic clusters was already present in the classical Schelling model and thus is mainly 
due to the behavioral rules of the latter. We can therefore imagine the process at play in our model as the combined 
outcome of a local birth and death process (migration in-out) on the single node and a coarse-grained node dynamics 
on the lattice. 


IV. AN ANALYTICAL SIMPLIFIED MODEL 

To gain insight into the behavior of the previously presented model, we hereby introduce a model ables to capture 
the main behavior of the metapopulation Schelling model, but simple enough to be analytically tractable. Because we 
consider the case for extremely large L, our model complements the one proposed in m devoted to the case L = 2. 

Using the notations introduced before the state of the system is thus completely characterized by the knowledge of 
(n^(t), n^(t)), being n^{t) = {nf{t ),... ,nfj{t)) and n^{t) = (nf (t),... ,nf ^). For a sake of simplicity we decided 
to compute the fitness using only the information from a single node, Eq. with j = i. The system evolution is 
done as previously and we still use the weak liquid version, an unhappy agent will move to an uniformly randomly 
chosen new node, provide there is enough space there. 

The model is thus intrinsically stochastic and hence it can be described by the probability P(n^,n^,t) to be at 
time t in state (n^,n^), whose evolution is dictated by the master equation: 

P(n^, n^, t + 1) = P(n^,n^, t) + (6) 

[T{n^,n^\n^ ,n^ )P{n^ ,n^ ,t) 

{fi^' ,n ^') 

— T{n^ n^)P(n"^,t)] , 


2 


Excluding the case where there is not strict majority in a node, this definition can be restated as: i and j belong to the same cluster if 
{nf -nf){nf - nf) > 0. 
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FIG. 2: Global behavior. Panels A: Nodes occupancy at equilibrium, each node has been colored according to the majority of 
agents in it: white empty nodes, blue only B agents, R only red agents and shades for intermediate cases. Panel B: Gluster size 
distributions for several values of e = 0.3, 0.4, 0.5, 0.6, 0.7. Panel G: First and second largest cluster as a function of e. Panel 
D: Interface size as a function of time for several values of e = 0.3, 0.4, 0.5, 0.6, 0.7. Results presented in panels B, G and D are 
based on 100 replicas of the simulations. 


being T{n^ , n^) the Transition probability to pass from state (n^, n^) to the new compatible one ). 

For incompatible states we set T(n^ = 0. 

The non zero transition probabilities can be computed using the rules previously defined (see Appendi:?jB] for a 
detailed account of the transition probabilities calculation), for instance the transition probability that an A agent 
moves from node i-th to node j-th. is given by: 

Ti{nf -l,nf + l\nf,nf) = 

1 nf / nf \ ^ ~ ~ 

N nf + nf ~ \nf + nf y L 

The function 0(x), defined to be 1 if x > 0 and 0 otherwise, translates the willingness of an unhappy agent to move; 
smoother versions can be used as well. Let us observe that to lighten the notation we wrote only the variables whose 
values change because of the transition. 

From the master equation one can compute the time evolution of some relevant quantities, for instance the average 
number of X-agents in every node i at time t, {nf){t) = t), X = A^B. 

Using the expression for the transition probabilities Eq. [B3] and assuming correlations can be neglected, i.e. 
{(ntf) ^ (n^)^, we obtain a system of finite differences describing the evolution of (nf) and (nf). 

To go one step further, we introduce the average fraction of A and B in each node, ai = {nf)lL and Pi = {nf) /L^ 
we rescale time by defining s = t/L and finally we assume each node to have an infinite large carrying capacity (see 
Appendi:^ for a more detailed discussion): 

daijs) ^ _ aj Q / Pi 
ds ai pi \ + pi 

1 — aj — Pi aj 
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dpi{s) 

ds 


-P 


Pi 

Oii + Pi 


0 




l- Uj- Pi pj 
N ^aj+P, 



(8) 


where p is the fraction empty nodes, (1 — ai{s) — Pi{s)) = pN. 

To disentangle the evolution of ai and Pi we introduce a new set of variables, the node emptiness^ 7 ^ = 1 — ai — Pi^ 
and the difference of fractions of A and 5, Q = — Pi. The new variables range in Q G [—1,1] and 7 i e [ 0 , 1 ]. 

Introducing the functions 


F{x) = x0(l — X — e) + (1 — x)Q{x — e) 
G{x) = xQ{l — X — e) — {1 — x)Q{x — e), 


we can rewrite Eqs. (B7) and iBSh as follows 


dlijs) 

ds 


pF 



2(1-70/ 



2(1-70/ 


and 


dCijs) 

ds 



2(1-70/ 




V 

2(1-7i)/ 


(9) 


( 10 ) 


The agents initialization used in the previous section, namely pN density of vacancies and an equal density {l—p)N/2 
of A and B agents, translate into {Ciili) uniformly distributed in a neighborhood of (0, p). We thus divide the domain 
of definition of into four zones (see Fig. and we will look closely to the dynamics in the Z 2 zone if e < 1/2 

and in the Z 4 zone if e > 1 / 2 . 




£> 1/2 


FIG. 3: The four zones Zi used to study the solutions of system (H, 


Let e < 1 / 2 , then one can easily prove that if (C, 7 ) € Z 2 one has r’ (5 + 2 ( 1 - 7 ) ) 
SO assuming that for all i one has (Ci(0),7^(0)) G Z 2 , then Eqs. (0 and ( p[Q| ) rewrite: 


= 1 and G 


(2 + 2 (iG)) 


djijs) 

ds 

dCih) 

ds 


= P-Ii 


-Pi 


Ci 

-li 


7i V 

N ^3 l-7j ’ 


C 

1—7 ’ 


( 11 ) 


as long as (Ci(t), 7 i(t)) will not leave Z 2 . Assume for a while this statement to hold, then the first equation can be 
straightforwardly solved to give 7 i(t) = p + e“^( 7 G 0 ) — p), that is for all i, 7 i(t) ^ p when t ^ 00 . The second 
equation can also be solved (see Eq. [S13]) and thus to prove that Q{t) 0 for all i when t ^ 0. This proves also a 
posteriori that {(i{t), ji{t)) will never leave Z 2 . 

The average magnetization rewrites in such variables as: 


(m) = 


7 w 

iv ^ 1 - 77 


( 12 ) 
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we have hence proved that for initial conditions in the magnetization asymptotically vanishes (see Fig. |^. 

The remaining case, e > 1/2, can be handle as well but it is more cumbersome (see Appendi:?|^- Let us only 

observe there that for (C^t) ^ ^4 one has F ^ 2 ( 1 - 7 ) ) ~ ^ ^ 2 ( 1 - 7 ) ) ~ assuming that for all i one has 

(Ci(0), 7i(0)) G Z4, then the right hand side of Eqs. © and ( [1Q| ) identically vanishes and thus the average magnetization 
depends on the domain where initial conditions have been set. However for a fixed size of the latter, one cannot satisfy 
the hypothesis (Ci(0),7^(0)) G Z 4 if e is closer enough but larger then 0.5, indeed the Z 4 zone shrinks to zero in this 
case (see Fig. |^. In this case one should take into account the dynamics of orbits whose initial conditions are set in 
Zi and Z 3 (see AppendiJSand Figs. 8 and 9). In conclusion the analytical model perfectly fits with the ABM for 
e > 0.5 while it has a different behaviour for e < 0.5, because in this range the ABM dynamics is strongly dictated 


by the stochasticity of the model, orbits tending to converge toward the equilibrium ( 0 ,p) ((/r) 
by fluctuations and thus sent into the zones Zi and Z 3 ((/i) ^ I). 


0 ) are destabilized 





FIG. 4: The average asymptotic magnetization as a function of e. Each point is the average over 50 simulations whose initial 
conditions are close to the equilibrium point — p and = 0. Black circles represent the results of the analytical model while 
grey diamonds to the ABM with 1-node fitness. Parameters are: A = 100 and p = 0.9. 


V. CONCLUSIONS 

We proposed and analyzed an extension of the classical Schelling model to a metapopulation framework, whose main 
outcome is the spontaneous emergence, for low values of the tolerance threshold, of heterogeneously populated nodes 
without any exogenous preferential attachment mechanism (for ID lattice this phenomenon induces the formation of 
urban skylines, Eig.|^. This behavior is connected to the permanence for long time of the system in a quasi-stationary 
non segregated state, where in each node the two populations are equally distributed. This quasi-stationary state 
can be recovered as stable equilibrium of the simplified analytical model. We can evince that the system stabilization 
toward the magnetized state (for the ABM) passes through the creation (by random agents moves) of highly populated 
nodes, hereby named towers. At the same time, global patterns emerge as in the classical Shelling model. Eigurej^ 
summarizes the possible behaviors of the system. Eor e < 0.5 the global clusters are described as neighborhoods 
formed by a strong majority of individuals of the same color (soft segregation). Eor e > 0.5 we have the formation of 
ghettos (hard segregation). 

The mechanism of formation of the clusters is a typical coarsening phenomenon for low values of e. On the contrary 
for low tolerance cases the clusters are formed around the towers that become stable points for a certain type of nodes 
(once an agent, whose kind corresponds to the majority already inside the tower, enters she never gets out). Once 
a higher density zone starts to exist, this mechanism reinforces the (majority color) population growth in this node 
and in the neighborhood. The global patterns start therefore to stabilize around the towers. 

Local magnetization phenomenon has been explained using an analytical approach able to describe the different 
equilibria of the model and the origin of the quasi-stationary state for low tolerances. 
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FIG. 5: Summary of the global outcomes of the Schelling metapopulation model on a ID lattice. Each box represents a node 
of the 1-dimensional lattice; the height of the box is given by the number of Blue and Red agent inside the node. The colored 
rectangle below the nodes represents the composition of the node obtained using the same color code used in Fig. Panels 
A1-A5. 
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Appendix A: Patterns dynamics 

The aim of this section is to present some details concerning the convergence of the agent based model to the 
asymptotic equilibrium pattern. We start by considering the convergence time, defined as the time needed for the 
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system to reach such equilibrium where no agents has incentive to move anymore, and its dependence of the tolerance 
threshold e. Intuitively, if e is large then most agents are happy and thus they will not move, hence the equilibrium will 
be reached quite soon. On the other hand if e is small, agents will be very often unhappy and thus they will relocate 
themselves to reduce this uneasiness, this will increase the convergence time. In the limit e ^ 0 this equilibrium will 
be never achieved and thus the convergence time will diverge. 



FIG. 6: Convergence time as a function of e. The blue line represents the average convergence time as a function of e for 100 
replicas of the system with same initial conditions and parameters values. The inset shows the boxplots for the convergence 
times for some representative values of e = 0.3, 0.4, 0.5, 0.6, 0.7. 

In Fig. we report the results obtained by the simulation of the metapopulation model and we can observe that the 
convergence time does not have a monotonic behavior as a function of e. Passing from from e = ltoe^0.6we observe 
an exponential growth of the convergence time (note the logarithmic vertical scale), these values of e correspond to 
the absence of macroscopic clusters (see Fig. 5 main text). Let us notice that the time to converge increase between 
e ^ 0.5 and e ^ 0.6 is due to a slower dynamics for the latter case (to move is less probable). Then for 0.3 < e < 0.6 
the convergence time exhibits an almost stable value - horizontal plateau - with a local minimum at e = 0.5, this is the 
phase where macroscopic clusters are formed through a coarsening process (see Fig. 5 main text). Finally for e < 0.3 
the convergence times start to increase faster than exponentially, the dynamics exhibits quasi-stationary states and 
the system reaches the equilibrium through the formation of towers (high densely populated set of contiguous nodes). 
Let us observe that asymptotic equilibrium is affected by the intrinsic stochasticity of the ABM, the convergence time 
will thus reflect this fact by showing a large variance, mainly for small e (see inset of Fig. [^. 

The patterns arising for small tolerance threshold are intriguing and peculiar to our model. As already observed 
the system spends quite a long time in a quasi-stationary (almost) homogeneous state and then suddenly jumps to 
a macroscopic segregated one. Such behavior is schematically represented in Fig. in the case of the ID-lattice 
for e = 0.25. The system stabilization toward the magnetized state passes through the creation (by random agents 
moves) of a highly populated nodes, hereby named towers corresponding to monochromatic clusters in the 2D model 
presented before. Such towers become attracting selective places for each kind of agents: once an agent, whose kind 
corresponds to the one of the majority already inside the tower, enters she never gets out from the tower because she 
will be happy there and this will increase the unhappiness of agents of the opposite kind still inside the tower. The 
net result is that once a highly dense zone starts to exist, this mechanism reinforces the (majority color) population 
growth in this node and in the neighborhood, because of the way the fitness is computed, using nodes at distance 1. 
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The global patterns start therefore to stabilize around the towers. As the first kind of agents stabilizes, with a certain 
delay also the stabilization of the second population is reached. Once the towers/clusters are formed, the complete 
equilibrium is reached as the interface between Blue and Red zones becomes empty. The lower is the tolerance value 
e, the higher is the initial density unbalance needed to start the stabilization process. 



FIG. 7: Generic time evolution toward a hard segregated pattern. Upper main plot: first and second cluster size as a function 
of time. Lower main plot: average magnetization of the system as a function of time. Left and right panels: System snapshots 
at different times. Each box represents a node of the 1-dimensional lattice, the height of the box is given by the total blue and 
red population of the node. The colored rectangle under the cell represents the majority color of the cell. 


Appendix B: More details about the analytical model 

The goal of this section is to present more details of the analytical model we introduced to capture the main behavior 
of the metapopulation Schelling model presented in the main text. 

Let us denote by respectively nf{t)^ the number of agents of type A, respectively 5, at node i-th at 

time t. Each node is also characterized by a number of vacancies nf{t) at time t, however for all t and all i one 
has nf{t) + nf{t) + nf{t) = L and moreover the number of agents of each kind is a preserved quantity, i.e. the 
system is closed. The state of the system is thus completely characterized by knowledge of (n^(t),n^(t)), being 
n^(t) = and n^{t) = {nf (t),... 

An agent A in node i-th is unhappy - or her fitness is low - if the fraction of B in the same node, that is nf / (n^+nf), 
is larger than a given tolerance threshold e G (0,1): 

A is unhappy at i if . ^ ^ > e , (Bl) 

nf + nf 


and similarly for B. 

We assume that once an agent decides to move because unhappy, she will move to an uniformly randomly chosen 
new node, provide there is enough space there (that we named weak liquid Schelling model). The model is intrinsically 
stochastic and hence it can be described by the probability to be at time t in state (n^,n^), that is P{n-^,n^A)‘ 
The evolution of such probability can be obtained using the master equation: 


P(n"^,t + 1) = P(n"^,t) + T(n"^,)P(n"^ ,,t) — n^)P(n"^,n^,t) (P2) 


{fi^' ) 
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being T{n^\n^'\n^^n^) the Transition probability to pass from state to the new compatible one, i.e. the 

system can pass from the former to the latter, Non compatible states cannot be linked and thus we will 

set T(n^ ,n^ = 0. 


n ,n |n ,n 

The non zero transition probabilities can be computed using the behavioral rules of the model: 


A moves from node i to node j: Ti{nf — 1, '^f) = 


1 


N nf 


r0 


nf + nf 


n 


E 


A moves from node j to node i: T 2 {nf + 1, — l|nf, nf) = 


B moves from node i to node j: Ts{nf — l^nf l\nf ^nf) = 


B moves from node j to node i: T 4 (nf + l,^f — l|nf ,n^) = 


N nf ^ nf 

1 nf 

N nf nf 

N nf + nj 


e 


e 


e 


A , R 

nf + nf 


A , R 

nf + nf 


nf+nf 


— e 


— e 


— € 


L 

E 


n 


\ ■ (B3) 


Let us observe that to lighten the notation we wrote only the variables whose values change because of the transition. 
The function 0(x) is defined to be 1 if x > 0 and 0 otherwise; smoother versions can be used as well. Being the 
structure of such transition probabilities very similar we will detail the way we compute the first one: 

• ^ : probability to draw the i-th node with uniform random probability; 


nf+nf 

. 0 ( 


: probability to draw one A agent among the nf present over the total node population nf + n 


B. 

i ’ 


nf+nf 


: probability the selected A agent is unhappy; 


• probability to select a vacancy in node j. 


The master equation is unmanageable but one can go one step further by computing the time evolution of relevant 
quantities, for instance the average number of agents A and B in every node i at time t, {nf) (t) = n^P(n^, t): 

{nf){t + l) - (nf)(t) = y;y; [nfTi(nyn/|nf+ l,n/-l)P(nf+ l,n/-l,i) 

n^ 

+ nfT2{nf,nf\nf — l,nf + l)P{nf — l,nf + l,t) — nfTi{nf — l,nf + l\nf ,nf)P{nf ,nf ,t) 

- nfT2{nf + l,nf - l\nf ,nf)P{nf ,nf ,t)\^ 

= -^{Tx{nf -l,nf + l\nf,nf)) + ^{T2{nf + l,n/ - l\nf,nf)). 

And similarly for B: 

{nf){t + l) - {nf){t) = -^{TAnf - l,nf + l\nf ,nf)) +^{TAnf + l,nf - l|nf,nf)). 


Using the expression for the transition probabilities (B3) and assuming correlations can be neglected, i.e. {{nf)‘^) 
{nf)‘^^ we get: 


{nf){t + l)-{nf){t) = / B\ ^ ( / / B\ 

TV ^ (nf) + (nf) V(^f) + {nf) J 


1V 

N 


(nf) 


-0 


(nf) 




and 




{nf) 


-0 


Nj+,{nt) + {nf) \{nt) + {nf) 


(nf) 


— e 




inf) 


-0 


«) 


, L-(nf) 

- {nf) 

/ L 


L-{nf)- 

{nf) 

L 


\ L - {nf) 

- {nf) 

/ L 


L - {nf) - 

{nf) 


(B4) 


^ f+i {^f) + {A^ V 


L 


(B5) 
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Observe that the right hand sides of (B4) and (B5) remain unchanged if we allow both sums to run over all node 
indexes, namely include also j = i. 


To go one step further, let us define the fraction of A and B in each node, that is 

ai = —^ and I3i = —^— , 


(B6) 


then we can rewrite ^ the previous equations (B4) and (B5) in terms of and Pi. Finally rescaling time by 5 = t/L, 

dividing the equations for ai and Pi by 1/L and passing to the limit L +oo we get: 


dai{s) _ ai{s + 1/L) - ai{s) 

ds L —>--|-cxD \ j L 


ai 


0 


N ai -\- Pi \ai P^^ 


Pi 


Ed 


ai 


' Pj) + 


ai 


■A 


N 


- E Q 

^ ai + /?,■ 




aj + pj 


(B7) 


and 


dpijs) 

ds 


lim 

L —>--|-CXD 

1 


Pi{s + l/L) - Pi{s) 


1/L 


A 


e 


N ai + Pi \ + /3; 


ai 


Ed 


a, i 


idj) + 


1- ai- Pi 

N 


V- 


A 


aj + P^ 


0 


ap 




(B8) 


Remark B.l (Some preserved quantities). Let us observe that the model eorreetly preserves the total fraetions of A 
and B agents in time and thus the total vaeaneies (1 — — Pi{s)) = (1 — (^^(O) — A(0)) = y where p is 

the emptiness defined previously, i.e. the total fraetion of vaeaneies. 

To pro ve this statement is enough to take the time derivative of ^^{1 — afis) — Pfis)) and observe that using 
Eqs. (B7) and (B8) one gets: 


— E(i -«id) - A(s)) = 0. 


One ean similarly prove the statement about the total fraetion of A and B agents 
So in conclusion the system is ruled by the following system of differential equations: 


daijs) _ 
ds t' 

dPiO) _ Pi 

ds cxi+Pi 


P ( hi — A jl-ai-Pi) ^ aj p f 

^ {ai+Pi 7 ^ N 2^j aj+pj ^ 

—Si_d I dj 0 / 

\ai+Pi N 


0 




— e 


aj+Pj 


(B9) 


Appendix C: The average magnetization 

To better analyze the system and in particular be able to describe the dependence of the average magnetization on 
e, we introduce a new set of variables, the loeal emptiness, ^i = 1 — ai — pi, and the loeal differenee of fraetions of A 
and B, (i = ai— Pi. The original coordinates can be obtained back using ai = {l—ji^(i)/2 and Pi = {l — ^i — fi)/2. 
The new variables ranges are Q G [—1,1] and 7^ G [0,1] and the magnetization rewrites as 


Let us introduce the functions 

F{x) := x0(l — X — e) + (1 — x)Q{x — e) and G{x) := xQ{l — x — e) — {1 — x)0(x — e), (C2) 


^ Let us observe that (nf) / {{nf) + {nf)) — e is positive if and only if OLifioLi + di) — e > 0, and similarly for the other term. 






































then observing that 
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OLi 


O-i + Pi 


Ci 


and 


A 


2(1 — 7i) + P'. 


Ci 


2 ( 1 -7i)’ 


we can rewrite Eq. (B 9 ) as follows 

d7i(£) ^ pp ^1 



f'C (3 + 5a5^) + F Ej C (3 + 3(1-1,)) ■ 


2(l-7i) 


(C3) 


To qualitatively study the solutions of the latter system we define four zones in the (C, 7) plane as follows (see Fig. 
3 main text): 


Z3 = {«.,)c|-l.llx|0.11)l-5ij^ 

Z 3 = {(C.t) e 1-1.1] X |0,1] : - - 2(I~) 

Z. = {«.7)U-l.l]x|0.1])l-j^ 


0 


2(1 - 7*) 


— e > 0 and - 


— e > 0 and ^ ^ —- 

2 2 ( 1 - 7 ,) 

1 Ci 

— e < 0 and —^,-r 

2 2 ( 1 - 7 ,) 


— e < 0 and - 
2 


Ci 


2 ( 1 -7z) 


-e < 0} 
-e > 0} 
-e > 0} 
-e < 0}. 


Let us observe that if e < 1/2 then Z4 is empty and if e > 1/2 then Z2 is empty. 

The initial condition of the Shelling model presented in the main te xt tr anslate into (Ci^Ti) distributed close to 
( 0 , p). We will thus be interested in the behavior of the solutions of Eq. (C 3 ) in the Z2 zone if e < 1/2 and in the Z4 
zone if e > 1/2. 

Let us consider the case e < 1/2, then one can easily prove that for all (C,7) ^ ^2 one has E + 2(1-7) ) ~ ^ 

G (1 + 2(1^7)) ~ assuming that for all i one has (C,(0), 7 ,( 0 )) G Z2, then: 


dCijs) _ Ci \ j±sr^ _k 
ds ~ ^l-7i ^ N 1- 


0 

J l-7i ’ 


(C 4 ) 


as long as (C,(t), 7,(t)) will not leave Z2. Assume for a while this statement, then the first equation can be straight¬ 
forwardly solved to give 7,(t) = p + e“^(7,(0) — p), that is for all i, 7,(t) ^ p when t ^ 00. 

The solution for the second is more cumbersome, but one prove the existence of a lower (upper) solution C^{t) 
(Ctd)) such that < Ci{t) < 7 (t) where: 


CHt) = 


e-*(l-7i(0)) 

1 -p-e-*(7i(0) -p) 


. 3X37 3. 37 3 n ^ fP+e-pim-p) x{l-x)^P ' 

Ci( 0 ) T 7 2e) ( j / „ dx 

\l-7,(0)y Jji(o) (x-p)i-p 


(C 5 ) 


from which one can get (i{t) 0 for all i when t ^ 0. 

Because the point ( 0 ,p) belongs to Z2 we have a posteriori proved that (C,(t), 7,(t)) will never leave Z2. A simpler 
but also weaker statement can be obtained by observing that the equilibrium = (0, p) for all i is stable being 

its eigenvalues —pi{I — p) and —1, each one with multiplicity thus there exists a neighborhood of (0,p) such that 
all orbits whose initial conditions are inside never leave it. Using the definition of (p) given by Eq. (Cl), we have thus 
proven that for e < 0.5 the average magnetization converges asymptotically to 0 (see Eig. 4 main text). 


Assuming now e > 1/2, then one can easily prove that for all (C, 7) G Z4 one has F 2{i-'y ) ) ~ ^ F 2(1^^^) ) “ 
0, so assuming that for all i one has (C,(0), 7 ,( 0 )) G Z4, then: 


’ dli{s) _ 
ds ~ ^ 

dCijs) _ PI 

d.'n ’ 


(C 6 ) 


and thus (C,(t),7,(t)) will not evolve and thus remain in Z4. However fixed the size of the domain centred at 
(C,7) = (0,p), where initial conditions are taken, the above assumption cannot be satisfied if e > 0.5 but close to it. 
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Indeed looking at the right panel of Fig. 3 in the main text, we realize that the zone Z4 shrinks to zero as e ^ 0.5 
from above. The aim of the rest of this section is to prove that we can take this fact into account and explain the 
behavior of (/i) also for e > 0.5. 

As already remarked the functions F and G vanish for (C^t) ^ ^ 4, thus Eq. (C3) can be rewritten as: 


Cj 


(2 + 2(l-7i)) w [SjeZs (2 + 2(l-7j) 


) + Ejezi ^{1 + 2(1 A)) 


ds 


+ N [SieZa (5 + 2(T^) + EjeZi G' (| + 2(1^) 


(C7) 


where we used the shortened notation j G to mean {Cj^lj) ^ and where the zones Zi and Z^ have been defined 
in the Fig. 3 of the main text. 

The initialization of the metapopulation model, translates into original variables initially uniformly randomly 

distributed in (1 — p)/2 + 1,1], for some positive small 5. Such domain is distorted passing to the new variables 

{zetai^^i)^ in particular it is translated into (0,p), rotated by 45° and expanded by a factor ^/2 in both directions 
(see left panel of Fig. 


Let us observe that F(0) = G(0) = F{ 1 ) = G(l) = 0 and thus Eqs. (C7) admit as equilibrium points Q = ±(1 — 7i). 
In conclusion, initial conditions inside Z 4 will remain there while orbits originated from initial conditions in Zi and 
Z3 will converge somewhere onto the straight lines C = ±(I — 7). In the left panel of Fig. we report 5000 initial 
conditions built using the above describ ed p rocedure for 6 = 0.02, points in Zi are marked in green, points in Z4 in 
blue and points in Z3 in red. The Eqs. (C7) are then numerically solved and the asymptotic positions are drawn in 
the right panel of Fig. giving to any points the color corresponding to the zone from where is started. One can 
clearly observe that all blue points remain in the Z4 zone, while green (respectively red) points converge to C = 7 — I 
(respectively to C = I — 7 ). 



FIG. 8 : Dynamics for e > 0.5. Left panel: 5000 initial conditions are uniformly drawn {ai,/3i) G (I — p )/2 + SU[—1,1] and 
then transformed into the new variables ( 0 , 7 ^); points in Zi are colored in green, points in Z 4 in blue and points in Z 3 in red. 
Right panel: asymptotic configuration of the orbits corresponding to the initial conditions given in the left panel, each point is 
drawn using its initial color. The black solid line is the straight line ( = (2e — I)(I — 7), the dot-dashed line is the straight line 
C = (I — 2e)(I — 7) and the dashed lines (right panel) are the straight lines ( = ±(I — 7). 


We are now able to provide a good approximation of the average magnetization. First of all let us rewrite Eq. (Cl) 


(/^) = 


I 

N 


5 : t7) 

deZ4 ieziuza 


^ ^ iCil -^13 1 ^ iCil 


where A^4 is the number of points in the Z4 zone and ^"13 the total number of points in the Zi and Z3 zones. Because 
points in Zi and Z3 converge to ( = — the right most term in the previous equation is trivially equal to A^i3/A^. 

To compute the contribution arising from points in Z 4 is more cumbersome, let us hereby stress the main ideas. One 
can easily show that if e > 0.5 +(5/(1 — p) then the lines ( = (2e — I)(I —7) and C = (1 “2e)(I — 7) do not intersect the 
diamond like domain (see left panel Fig. and thus all points are initially in the Z 4 zones, this means that N 4 = N 
and A^i 3 = 0 so the magnetization has only a contribution from the Z4 zone. On the contrary if e < 0.5 + s/{l-p) 
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e>l/2 + 8/(l-p) 


e<l/2 + 8/(l-p) 


FIG. 9: The geometries used to compute an approximation to the average magnetization in the case e > 0.5. 


some initial conditions are in the and Z3 zones (see right panel Fig. and thus the magnetization has both 
contributions. 

In the case e > 0.5 5/(1 — p) one can estimate the contribution to the magnetization as: 




^ —C+p+2(5 


C+p—2(5 




1 

^Jo ' 


rfCCiog 


1 - c - P + 2(5 

1 + C — p — 2(5 ’ 


(C8) 


being 8(5^ the measure of the blue diamond domain in the left panel of Fig. While for e < 0.5 -1- 5/(1 — p) one can 
find the following estimation 


= 2_ W r^—[ 


^(Cp+Cq)/2 c+p+2(5 -I or 

dCC dj- -= -/ 

Jc+p-^d 1 — 7 V Jo 


(Cp+Cq)/2 1 _ ^ ^ I 2(5 

«=») 


where V is the measure of the polygon TPQRQ'P' (in blue in the right hand side of Fig. |^. Let us observe that 
in the last integral we make the approximation that points P and Q can be replaced by an average point whose ( 
coordinate is the average of the ones for P and Q; this is a minor assumption that helps to compute the integral and 
will not influence the final result as show below. Now both integrals can be exactly computed. 

To get the final estimate for (p) we need to compute N4/N and Ni^/N in the case e < 0.5 + (5/(1 — p), the other 
case being trivial and already considered. Assuming the number of points sufficiently large we can affirm that such 
fractions are well approximated by the ratio of the corresponding polygons, that is 


A^4 V 

N 8(52 


and 


N ^ 8(52 ' 


In conclusion we can approximate the average magnetization for all e > 0.5 with the formula: 




(CIO) 


where pint 1 is given by Eq. (C8) valid for e > 0.5 +J/(l — p), and Pint,2 by Eq. (C9) for the case e < 0.5+ (5/(1 — p). 
In Fig. 10 we compare th e av erage magnetization obtained using its ver y firs t definition Eq. (Cl) and the numerical 

with the approximation given by Eq. (CIO) for several values of 5 and p = 0.9. 


integration of the system 

One can observe the very good agreement, providing thus an a posteriori validation of the assumptions made so far. 
Comparing with Eig. 4 of the main text, this provide also a very good approximation to the average magnetization 
computed using the metapopulation model. 
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FIG. 10: Averaged magnetization {/j ) an d the approximation ^ function of e > 0.5. Symbols correspond to 

numerical integrati ons of the system (C3) and the use of the definition for the magnetization, lines to the calculation of the 
approximation Eq. (CIO). Blue circles correspond to (5 = 0.02, black squares to (5 = 0.01 and red diamonds to (5 = 0.005. Each 
point correspond to the average over 50 replicas. The emptiness has been fixed to p = 0.9. 












